Fluid-solid transition in hard hyper-sphere systems. 
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In this work we present a numerical study, based on molecular dynamics simulations, to estimate 
the freezing point of hard spheres and hypersphere systems in dimension D = 4, 5, 6 and 7. We 
have studied the changes of the Radial Distribution Function (RDF) as a function of density in 
the coexistence region. We started our simulations from crystalline states with densities above the 
melting point, and moved down to densities in the liquid state below the freezing point. For all the 
examined dimensions (including D = 3) it was observed that the height of the first minimum of 
the RDF changes in an almost continuous way around the freezing density and resembles a second 
order phase transition. With these results we propose a numerical method to estimate the freezing 
point as a function of the dimension D using numerical fits and semiempirical approaches. We 
find that the estimated values of the freezing point are very close to previously reported values from 
simulations and theoretical approaches up to D = 6 reinforcing the validity of the proposed method. 
This was also applied to numerical simulations for D = 7 giving new estimations of the freezing 
point for this dimensionality. 
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I. INTRODUCTION 

Although the study of fluids of hard spheres in dimen- 
sions higher than three is not new-i~— it has attracted re- 
newed attention in recent yearsi This interest is because 
as we change the dimension these systems share common 
properties that may lead to general relations among them 
and eventually could bring new insights to traditionally 
difficult problems for two and three dimensional systems. 

One of such properties is the fluid-solid transition, 
which has been widely studied in hard-sphere fluids and 
is known to appear also for hard-hypersphere fluids. Ev- 
idence of this transition found by numerical simulations 3 
and theoretical estimations 8 suggests that it preserves 
the same characteristics of the three dimensional case. 
Nevertheless, the exact determination of freezing and 
melting points at any dimension (including two and 
three) is, at present, not possible. Only using numeri- 
cal simulations and theoretical approaches is how these 
values are known with some accuracy. 

The freezing point of the hard-sphere (HS) fluid was 
first observed and estimated by Alder and Wainright 
in 1957 using molecular dynamic simulations^—. Since 
then the problem has been revised using different simula- 
tion techniques, ranging from Monte Carlo methods 1 ^ - — 
to the recent Direct Coexistence Simulations 14 . Also 
different theoretical approaches have been applied to 
the same problem, such as Mean-Field Cage Theories 
(MFCT) 1 ^ and Density Functional Theories (DFT)i£~— , 
in particular the fundamental measure theory (FMT) 
by Rosenfeld^ describes both the fluid and solid phases 
while keeping the correct limit in the close packing. An- 
other advantage is that, in contrast to other density func- 
tional approaches, it does not require the direct correla- 
tion function as an input. All of the above works seem 
to coincide that freezing occurs at a packing fraction r\ 
close to t]f ~ 0.494. This observation has been also 



confirmed by different experiments with colloidal suspen- 
sions of hard-core particle a 21 i 22 . 

For hard-hypersphere(HHS) fluids at dimensions D > 
3, the information available is so far limited. In 
1984, Leutheusser 2 found an analytic solution, using the 
Percus-Yevick approximation for the HHS problem that 
allows one to obtain algebraic expressions for the pres- 
sure equation of state in odd dimensions. At present such 
equation of state has been computed only up to D = 7 
using the compressibility and the virial routesSS,. Unfor- 
tunately, the Percus-Yevick approximation gives no in- 
formation on the freezing transition and neither do other 
empirical and semiempirical approximations to the equa- 
tion of state. Recently, some theoretical approaches have 
been developed for the prediction of the freezing tran- 
sition as a function of D. For example, Finken et al£ 
generalized the scaled-particle theory to arbitrary D for 
the fluid phase and a cell theory for the crystalline phase 
to predict a first order freezing transition up to D = 50. 
In the same direction Wang 1 -, based on a MFCT, has 
obtained a simple expression for the freezing point as a 
function of D. At the moment, a full validation of both 
predictions based on computational simulations can not 
be done because of the lack of results. As far as we know, 
numerical simulations have been carried out to obtain 
information on the fluid equation of state for D = 4 to 
D = 93j23-.26 ; but t ne freezing points have only been esti- 
mated for D = 4, 52i2L2£ and recently Van Meel et alM 
provided more accurate estimations up to D = 6, while 
for D = 7—, the phase transition has been reported but 
its location has not been estimated. 

The main purpose of this work is to propose a simple 
empirical method, based on molecular dynamics simula- 
tions, to estimate quantitatively the freezing and melt- 
ing point for HHS fluids. Our proposal focuses on the 
changes in the Radial Distribution Function (RDF) when 
a fluid-solid transition takes place. For D = 2 and D = 3 
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it is well established that the RDF must contains relevant 
information on the structural precursors of freezing 3 ^ and 
therefore an empirical parameter defined on the basis of 
the changes in shape could be used to find the freezing 
point with some accuracy. Some empirical rules have 
been previously proposed for three dimensional fluids, 
like the Hansen- Verlet crystallization rule 3 -^, based on the 
height of the fist peak of the static structure factor. Also 
in 1978 Wendt and Abraham™, used the changes in the 
structure of the second peak in the RDF for metastable 
states to define a parameter that could allow them to 
estimate the glass transition density. Actually the prop- 
erties of the second peak of the RDF as density increases 
are directly related with the crystallization processes and 
the appearance of a shoulder on it is considered as a sig- 
nature of the fluid solid transition. The study of such 
changes for fluids at higher dimensions may give some 
criteria to estimate the freezing point on them, never- 
theless as dimensionality increases the oscillation of g(r) 
gets significantly dampe d 27 i 33 , and hence the character- 
istic shoulder of the second peak becomes difficult to ap- 
preciate. Based on this idea, we decided to use the prop- 
erties of the first minimum of the RDF instead, to find 
the freezing point for D > 3. 

The work is organized as follows: in section [XT] a re- 
view of the fluid-solid transition is included, as well as 
new molecular dynamics numerical simulations for the 
HS system at packing fractions close to and between the 
freezing and the melting points. This allows us to de- 
fine and validate the method we will later generalize to 
find the freezing and melting points of HHS systems for 
D > 3. In section Unl we report the results obtained with 
that numerical proposal for dimensions D = 4, .., 7. A 
discussion of the obtained results is included in section 
IIVI and the paper is closed with a summary of the results 
and further conclusions in section IVl 

II. FLUID-SOLID TRANSITION IN THE HARD 
SPHERE SYSTEM. 

The HS system is defined as a set of N identical 
molecules interacting with the hard-core potential: 

^(^) = {o %To 0) 

where a is the diameter of the spheres. 

Is well known that, because of the form of the in- 
teraction, the compressibility factor of the HS fluid: 
Z = PV/NKbT (where P is the pressure, N the num- 
ber of particles, V the system volume, T the tempera- 
ture and Kb the Boltzmann constant) can be written as 
a function of only one variable, either the reduced den- 
sity p = N<t 3 /V or the packing fraction r\ = pV sp h, with 
V sp h the volume of a unit diameter sphere. For a three 
dimensional HS fluid the volume of the unitary sphere is 
simply given by V sph =tt/6. 



The phase diagram can be reduced to only two sta- 
ble phases: fluid and crystal, defined completely by the 
freezing and melting packing fractions, rjp and t\m re- 
spectively. With some accuracy is possible to say that 
these packing fractions are located at rjp = 0.494 and 
rjM = 0.545 (see Table[T]) respectively. 

The pressure equation of state as a function of r\ 
draws a stable fluid branch and extends up to the freez- 
ing point, where the freezing pressure remains constant 
up to the melting point. At this point, the crystal 
branch appears extending up to the close packing den- 
sity (rjcp = 7r/3\/2). Compressing the hard sphere fluid 
beyond the freezing point but preventing crystallization, 
the system may enter into a metastable fluid branch. It 
has also been widely studied and it is known to drive the 
fluid phase to a supercooled liquid like states and undergo 
a glass transition. 

The freezing transition implies qualitative and quanti- 
tative changes in the structural properties, that may be 
observed in the Radial Distribution Function (RDF) . One 
characteristic effect observed in the HS fluid, is the ap- 
pearance of a shoulder in the second peak for dense fluid 
states close to the freezing point, which is associated with 
the formation of local crystalline regions 3 . . Another, per- 
haps less examined feature, is a change in the width and 
depth of the first minimum, that could be related with 
the cage effect produced by the local crystalline order. 

Our proposal here, is to examine using Molecular Dy- 
namics simulations the changes in the first minimum of 
the RDF with the purpose of defining a method to mea- 
sure, with reasonable accuracy, the freezing point in the 
three dimensional HS fluid. We believe that numerically 
this procedure may be simpler and more accurate than 
a method based on the evolution of the shoulder of the 
second peak. In addition this could also be simple to 
extended to hard hyper-sphere fluids at arbitrary dimen- 
sions. 



A. Simulation Details. 

In order to design a general simulation procedure com- 
mon to HS and HHS, we decided to start examining the 
three dimensional system, simulating states in the crystal 
branch from dense crystal to dense fluid phase. The posi- 
tions of spheres were set initially in a face centered cubic 
(FCC) structure with periodic boundary conditions and 
initial velocities chosen from a Maxwell Boltzmann dis- 
tribution. To establish a balance among simulation time 
and reliability of averages, we set the initial positions in a 
FCC crystal of 6 x 6 x 6 unitary cells, giving a total of 846 
particles. Since the number of particles and simulation 
volume cell are fixed, the diameter was used to change 
the packing fraction and therefore the density. In addi- 
tion, in this particular case, for comparison purposes we 
have used previous simulations for the metastable fluid 
branch which use 3951 particles in dense random packed 
states^. 
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Based on the algorithm of Alder and Wainwright^ 
to simulate the Hard Spheres fluid we performed event- 
driven molecular dynamics simulations, moving the sys- 
tem in successive collisions. The equation of state was 
measured computing the compressibility factor through 
the mean virial forces, given by the change on momentum 
at collision, as: 
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which is related to the compressibility factor through Z = 
1 + S. In Eq. ([2]) N is the number of particles, (v 2 ) is 
the mean square velocity, t is the time, N c (t) is the total 
number of collisions at t, f^j is the vector joining the 
center of the colliding particles i,j and Ai7, is the change 
in velocity of particle i. For the structural properties, 
RDF was obtained through the standard procedure by 
determining the mean number of particles at a distance 
between r and r + dr away from one of them, averaging 
over all particles and events. 

In order to obtain equilibrium properties, before the fi- 
nal simulation a stabilization stage was performed, where 
the system is left to reach the equilibrium ageing by 
one hundred thousand collisions and verifying that the 
change in the mean virial forces reaches a stable mean 
value. Once the system is stable the final configuration 
was taken as the initial condition for the true simulation. 
The final simulation stage was proposed to take averages 
of 100,000 collision steps. 



B. Equation of state, structural properties and the 
phase transition. 

The simulation results for the three dimensional case 
were validated using some analytic approaches for each 
branch. For the states in the fluid and metastable fluid 
phases we compared with the Carnahan-Starling equa- 
tion of stated. The crystal phase was approached with 
the Hall equation 3 -^ and for the glassy state we used 
the free volume approach from Speedy 3 -^. As can be 
seen in Fig. [T] numerical results are in good agreement 
with analytic equations, giving enough evidence that the 
simulation program is working properly. Although the 
used equations correspond to simple empirical expres- 
sions, they provide a good approximation to the general 
behavior of the system. Recent efforts by Bannerman 
et al.— have been made in order to obtain a new and 
more accurate equation of state, also compatible with the 
known high order virial coefficients. This reference also 
provides recent accurate compressibility data from simu- 
lations with large systems of about 10 5 and 10 6 particles, 
which are compatible with our simulation data (see Fig. 

It is worth to mention that, from the same figure, it is 
possible to note that molecular dynamics drive the sys- 
tem from stable crystal to metastable crystal states for 
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FIG. 1. Compressibility factor Z for the Hard Sphere sys- 
tem. The chosen simulation states for each phase are de- 
noted by circles (fluid), squares (metastable states) and di- 
amonds (crystal), simulation results by Bannerman (a) and 
the lines denote the equations of state: Carnahan-Starling 
(dot-dashed)^ , Speedy (dashed)^ and Hall (continuous)^. 



densities below the melting point. Such states are repro- 
ducible under the same simulation conditions, and the 
error bars which lie within 1% were not included to avoid 
overcrowding. 

Also for validation purposes the RDF was measured 
directly from simulations and compared with predic- 
tions from the Rational Function Approximation (RFA) 
method proposed by Yuste and Santos 39 . To com- 
pare with our simulations we used the Carnahan-Starling 
equation of state as input in the RFA scheme, predicting 
with excellent agreement the RDF's in the fluid phase 
but loosing accuracy close to the freezing point. 

Figure [2] shows some results for the RDF in the fluid, 
metastable and solid phases. For states just below the 
freezing point (77 = 0.47), the simulated RDF shows the 
characteristic shoulder on the second peak, indicating the 
rising of preceding structures of crystallization. Finally 
in the crystalline phase the RDF shows the rise of the 
peaks corresponding to the FCC structure, in such a way 
the simulations are giving expected results which can be 
used for more detailed analysis. 



C. The RDF analysis. 

The shape of the second peak of the RDF close to 
the freezing transition has been used in the prediction or 
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FIG. 2. Simulation (points) and theoretical (line) Radial distribution Function for densities corresponding to fluid, solid and 
metastable states. 



analysis of fluid-solid transitions—^. Nevertheless, the 
changes of the first minimum in the freezing transition, 
as far as we know, has not been deeply examined and can 
also carry important information. 

To measure the position of the first minimum from our 
simulation results, which we defined located at position 
r = X m i n and with depth g(X m i n ), we decided to fit the 
first minimum of the RDF with a polynomial function 
f(x) of the form: 

f(x) = — H h a 2 + a 3 x + a^x + a 5 x , (3) 

x A x 

to determine the coefficients aj by least squares and to ob- 
tain the minimum by solving the equation f'(X m in) = 0. 
The mean values of X m i n and g(X m i n ) were associated 
with the values obtained from for the final averaged RDF. 
The error bars were obtained computing the mean square 
difference between the instantaneous (at a time step, av- 
eraging only over particles) and the mean values in a 
sample of the last 1000 configurational states. 

Recently, using the short and long range asymp- 
totic behavior from the Percus-Yevick approximation, 
Trokhymchuk et al4^ derived an accurate analytical 
equation for g(r) 7 and the parametrized relations 

X m in = 2.0116 - 1.0647J7 + 0.0538?7 2 (4) 

and 

g{X min ) = 1.0286 - 0.6095/7 + 3.5781?y 2 (5) 
-21.365177 3 + 42.6344?7 4 - 33.8485r7 5 

for the location of the minimum in the range 0.1 < rj < 
0.47. These expressions were derived only for D = 3. 



The simulation results for X m i n and g{X m i n ) are 
presented as functions of the packing fraction in Fig. [3] 
and Fig. [4] For comparison, we have included the predic- 
tions derived from the RFA method and Trokhymchuk's 
expressions represented with continuous and segmented 
lines in both figures. Although the RFA predictions for 
X m i n are reasonably accurate while expression Q shows 
a better agreement, the result for g(X m in), is remarkable 
because the simulation data coincides with both the 
RFA theory and expression ([5]) for the whole fluid 
branch from low densities to the freezing point. Beyond 
the freezing point, the simulation data split into two, 
showing a clear difference between the metastable fluid 
and metastable crystal branches resembling a second 
order phase transition. It is important to remark that 
this change in g(X m in) may not be strictly interpreted 
as signaling a second order phase transition in a ther- 
modynamic sense i.e: equality of pressure and chemical 
potential in both phases. However, since for hard core 
systems the phase transition is purely entropic, one 
would expect that the changes in structure such as the 
one present in g(X m in) may indeed have a bearing on it. 
Therefore we assume that the changes on the mean value 
of Xmin and in particular g{X m in) are good indicators 
for the freezing transition and that the freezing point 
can be estimated using interpolation techniques on 
any of them. However, the use of g{X m i n ) looks better 
for that purpose, because the data show a cleaner aspect. 

We also studied the influence of the system's size by 
the inspection of the changes on measurements of X m m 
and g(Xmin) for a 2048 particles with respect to the 864 
particles system. Since the results obtained with 864 par- 
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FIG. 3. Comparison for the position of first mini muni ^Cmin of 
RDF for the hard sphere systems. The continuous line is the 
result from RFA method and the dashed line is the Trokhym- 
chuk et al. expression given in equation (Q. Diamonds (0) 
represent our simulation results with 864 particles and in the 
inset stars (*) correspond to 2048 particles. 
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FIG. 4. Behavior of the value g(X m i n ) with packing fraction 
for hard spheres: From the RFA method (continuous line), 
Trokhymchuk et al. expression given in equation J5| (dashed 
line) and simulation data for 864 (0) and 2048 particles (*). 



TABLE I. Estimated values of the freezing t\f and melting 
t/m packing fractions as well as the coexistence pressure for 
the HS system. 
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0.545 


11.7 
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0.487 


0.545 


11.54 
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0.491 


0.544 


11.54 
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0.476 


0.542 


10.1 


GELA 19 


0.495 


0.545 


11.9 


EXP.^ 


0.494 


0.545 




This Work 


0.488 


0.545 


10.95 



tides seems to reproduce well the theoretical estimations 
in the fluid branch, we choose to increase the system 
size for selected densities in the metastable crystalline 
branch, which is as well the most important for our pur- 
poses. The comparison is shown in the insets of Fig. [3] 
and Fig. 21 there we observed a maximal deviation of 2% 
on the mean value and a decrease in the error bars as we 
increase the number of particles more than double. Then 
we conclude that the shape obtained for both parameters 
will remain almost the same with an important raise in 
particle number. 

To estimate the freezing point from the simulation 
data, as a first approximation, we fit the simulation data 
of g(X m i n ) with two second order polynomials for the 
fluid and the metastable solid branches. The intersec- 
tion between these polynomials gives an estimation of 
the freezing point (tjf)- From the equations of state for 
fluid states, the corresponding coexistence reduced pres- 
sure is computed as p* F = rjpZ^p) /V sp h and the extrap- 
olation to the crystalline states determines the melting 
density t]m- These results were compared with previ- 
ously reported values in Table U showing that although 
the estimation is rough the technique gives comparable 
results. 

From this analysis we conclude that: 

1. Our simulations reproduce the stable branches 
of the phase diagram, as well as the metastable 
branches. 

2. The RFA predictions for the position of the first 
minimum give, compared with simulations, good 
qualitative and quantitative results for the fluid 
phase up to the freezing point. Beyond this point, 
it is not compatible with metastable branches. 

3. The analysis of the behavior of the first minimum 
on the RDF and the reasonable agreement of the 
estimated values with previous known results, sug- 
gest that these properties can be used to estimate 
the freezing point in an empirical way. 

In the next section we will discuss an extension of the 
numerical method, to approximate the freezing point for 
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TABLE II. Number of particles in a D-type lattice as function 
of the dimensionality D and the hypercubic cell size N c . 



TABLE III. Number of hypercells for simulation 



D 
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N c 
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32 


108 


256 
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128 


648 


2048 
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16 


512 


3888 


16384 
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32 


2048 


23328 


131072 


7 


64 


8192 


139968 


1048576 



systems at higher dimensions, examining the first mini- 
mum in the RDF of HHS from 4 to 7 dimensions, focusing 
on a quantitative analysis for the value g(X m i n ). 



III. THE EXTENSION TO THE CASE OF D 
DIMENSIONAL HARD HYPER-SPHERE. 

Based on the computational algorithm to simulate 
the hard sphere system^, we made an extension to 
D-dimensional hard hyper-sphere systems following 
the same simulation technique used by Michels and 
Trappeniers^ for D — 4 and D = 5. The computational 
code was designed to be extended to any dimension just 
by changing a library which includes all mathematical 
operations that depend on dimension. In this way the 
program use exactly the same dynamics on multidimen- 
sional spaces. 

For an arbitrary dimension D the generalization of the 
reduced density of a collection of hyper-spheres of diam- 
eter <7 is defined as p = Na D /V and allows us to define a 
generalized packing fraction by r\ = pV sp hD, where V sp hD 
is the volume of a hypersphere of unit diameter, which is 
given by: 



V sphD {a = 1) 



(tt/4) 



D/2 



n 



r(i + £>/2) 



where T denote the usual Gamma function. 

The simulation box used for all numerical experiments 
is an hypercube of side L with periodic boundary condi- 
tions and volume V = L D . To start from a crystalline 
state we choose to divide the simulation hypervolume 
in iV c D— type primitive hypercells of side I in each di- 
rection, such that V — (N c l) D . For a D— type lattice, 
the primitive hypercell contains No — 2 D ~ 1 spheres and 
therefore the number of particles in the simulation box 
as a function of N c and D is N p = 2 D - 1 N®. In Tabled 
we present the explicit numbers of particles in the simu- 
lation box for dimensionality between three an seven and 
N c = 1 to iV c = 4. Under this scheme the growth in 
the number of particles is geometrical and do the simula- 
tion very costly as the primitive cells and dimensionality 
increase. 

For the present simulations, and again to keep a rea- 
sonable balance between computing time and being able 
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32 
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32 
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to examine a wide density range, the initial configura- 
tions were chosen as shown in Table IIIII 

To validate our numerical results we considered recent 
simulation data from Skoge et alj^l, van Meel et al^S 
and Lue et al4i. The simulations of Skoge et al. were 
done for 4 and 5 dimensions with a maximum number of 
particles of 32768 and 124416, respectively. On the other 
hand the simulations of van Meel et al. were for 4, 5 and 
6 dimensions, using a maximum number of particles of 
4096 , 3888 and 2048 for the same D-lattices as initial 
states. Finally Lue et al. examined the solid lattices (-D4 
and D5) from event-driven molecular dynamics (10000 
and 16807 particles, respectively), and the metastable 
fluid phases by Monte Carlo simulations (from 4096 to 
10000 and 3125 to 7776 particles, respectively). Com- 
pared with these simulations the set of particles we used 
is very small. However it has been noted^ that, in the 
case of the compressibility factor at dimensions as high as 
D = 7, the comparison of numerical results with 64 par- 
ticles are in good agreement with simulations of systems 
of the order of thousands of particles. This may be con- 
firmed in Fig. [5j which shows the comparison with our 
simulation data. Size effects are very small for D = 4 
and D = 5. For D = 6 these changes seems to be small 
but may be important. Therefore, we also made a sim- 
ple examination of how the size of the sample affects the 
measurement of the value g(X m i n ) for the case D = 6, 
for which we use the smaller samples. Considering not 
only the hypercubic primitive cell (N p = 32), but also 
the rectangular simulation cells where one of the sides 
of the simulation cell is composed by 2 primitive cells 
(N p = 64) and the case where two sides have 2 primitive 
cells (N p = 128), to estimate the changes for a density 
close to the freezing point (p = 1.3). On the RDF one can 
notice first that, the larger sample is reflected in general 
as a smoother distribution. On the other hand, the evo- 
lution of the estimations of g(X m i n ) in time shows that 
the mean value converges with differences of the order of 
2% between systems of N p = 32 and N p — 64 particles, 



while the differences between N„ 



64 and N p = 128 are 



less than 0.5%. It seems that as we increase dimension- 
ality the growth of the system size affect less the equilib- 
rium properties. For D = 7 the comparison is possible 
only with simulation data in the fluid branch, however 
we believe the behavior of g(X m i n ) for larger simulation 
system will be close to the one we have obtained. 

For the whole range of densities, our values of the com- 
pressibility factor displayed in Fig. [BJ show for all ex- 
amined dimensions a similar qualitative behavior to the 
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FIG. 5. Comparison of measured Z (*) for set of particles 
given in Table Hill with data from Van Meel et al. (o), Skoge 
et al. (continuous line), Bishop et al. (□), Lue et al. (V) and 
Pades Zu^i (segmented line). 
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FIG. 6. Compressibility factor for D—4, 5, 6 and 7 as ob- 
tained by molecular dynamics simulations(dots) compared 
with reported semiempirical equations of state by Bishop et 
al.^^(lines) and fits (Eq. © and coefficients of Table HVJ) 
to the respective crystalline branches. 



D = 3 case. Simulation data agree with reported equa- 
tions of state up to D — 7 for the fluid phase. In par- 
ticular for Fig. [5] the continuous lines that reproduce the 
fluid phase, were chosen to be a Pade[4,5] of the form: 



1 



PiP 



P2P 2 + P3P 3 



l + qip + q 2 p 2 + q 3 p i + qAP + q^P 



(6) 



which is known to be accurate up to D = 9 using the co- 
efficients {pi} and {<7i} determined by Bishop et ah 25 ' 42 . 
For the crystalline branch, and for dimensions D > 3, 



TABLE IV. Fitting parameters for the equation of state Zc 
for the crystalline branch. 





Vcp 


PCP 


CO 


C2 


C2 


4 


0.61685 


2.0 


14.6452 


-18.6503 


9.07211 


5 


0.46526 


2.82843 


16.0717 


-16.0419 


6.36201 


6 


0.37295 


4.6188 


11.7468 


-7.02387 


2.46827 


7 


0.29530 


8.0 


17.3448 


-9.86783 


2.56618 



we propose that one accurate fit is given by: 



Zc(p) = 



Cp + cip + c 2 p 
1 - p/pcp 



(7) 



where pep is the close packed density. To be compati- 
ble with reported equations of state, this proposal is ex- 
pressed in terms of reduced density but may be rescaled 
to packing fraction. The fit of our simulation data gives 
the values of cq , ci and c 2 shown in Table IIV1 and the 
values of pep (or rjcp) were obtained from the lattice 
catalog of Nebe and Sloane 4 ^. 

Other works use the functional form by Speedy^ 4 - 



D 



cs P 



a ep (p/pcp - b sp ) 



1 - p/p 



CP 



P/P 



CP 



sp 



where the values of as P , bs p , and cs P are fitting constants. 
Both expressions are equally valid, and fit very well the 
crystalline branch. However, the proposed equation ([7]) 
seems to yield to better agreement for the metastable 
crystal branch. 

The structural properties are analog to the D = 3 case, 
showing the same disorder-to-order transition once the 
density increases and the system changes from fluid to 
metastable solid phase. 

In the analysis of the values of g(X m i n ) for D = 4 to 
7 shown in Fig. [3 it is remarkable that the analogy to a 
second order phase transition seems to be still valid for 
all examined dimensions. Therefore we applied the same 
numerical method used for D = 3 in the estimation of the 
freezing points for D — 4 to 7, i.e. through a numerical 
fit of two polynomials of order three to the mean values. 
The uncertainty was estimated taking into account the 
error bars of g(X m i n ). The results of this procedure are 
in Table PVTl under the label PF. The estimated values up 
to D = 5 could be compared with previous estimations 
by Michels and Trappenier 3 , Luban et al. 45 and recent 
results by Van Meel et ali22(see columns 1,2 and 3 in 
Table IVI[) . For all treated dimensions we compared also 
with the theoretical result of Xian-Zhi 15 , which are based 
on a mean field cage theory, and may be summarized in 
the general expression: 



PS 



p C p(l + 2D) 
2 + 2Dp CP a D ' 



(8) 



Again the values of pc p were taken from reference 4 ^. The 
agreement with our results is good until D = 6 and lie in 
the same order of magnitude up to D = 7. 
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IV. DISCUSSION. 

Following the work of Rohrman and Santos 46 , it is pos- 
sible to approximate the RDF of HHS fluids in odd di- 
mensions using a RFA method. Therefore by this method 
the changes on the first minimum of the RDF can be 
tracked as a function of the packing fraction and the di- 
mensionality. Unfortunately, this procedure is possible 
only numerically, but can give more information on what 
is observed in simulations. A numerical inspection of the 
predictions of the RFA method for D = 5 to D = 9 
suggests that the value of g(X min ) may be expressed by 
a function of rj and D which should comply with the 
following characteristics: i) in the limit of low packing 
fractions the value of g(X m i n ) must tend to 1. ii) The 
value of g(X m i n ) should decrease with density following 
a power law and Hi) the dimension modulates the width 
of the curve. A simple expression that meets these re- 
quirements is: 

g(X min ) = g M inL(v, D) = l- Drf, (9) 

where the coefficient a may depend on D. Using this 
equation to fit the numerical predictions of the RFA 
method for D=3,5,7 and 9 it is possible to establish that 
a must have a tendency that approximately follows an 
equation of the form: 

a(D) = a Q D- a \ 

where the two parameters involved best fit the predic- 
tions with the values ao = 1.59243 and ct\ = 0.5686. The 
advantage of deriving these relations is that expression 
^ may be assumed to be valid also for even dimensions, 
where no predictions are available from the RFA method 
and therefore may be used to compare all the simula- 
tion data obtained. Such comparison is shown in Fig. 
where the Eq. © was used to evaluate the continuous 
lines approaching the simulation points in the fluid phase. 
Clearly, they provide us the good estimations both qual- 
itatively and quantitatively for all of them. 

The value of g(X m i n ) in the crystalline branch can 
not be evaluated yet from any theoretical scheme, but 
examining with some detail the simulation data at even 
and odd dimension, one can note that the curves present 
an exponential decay with density. It is also reasonable to 
expect that g(X m i n ) goes to zero at some finite packing 
fraction lower than the close packing r\cp- Therefore, we 
propose the following general expression: 

9Minc(ji,D) =D- b ^+-)+d (10) 

where the parameters b, c and d depend on D and may 
be obtained as the best fit to the simulation data, giving 
the values included in Table IVl 

The changes observed for the parameters can be repro- 
duced with exponential functions of D. The best fitting 



TABLE V. Fitting coefficients for gMinc(T},D). 



d b c d 

3 7.559 -0.487 -0.425 

4 10.422 -0.296 -0.131 

5 20.715 -0.182 0.048 

6 30.999 -0.112 -0.049 

7 42.107 -0.067 -0.053 



TABLE VI. Comparison of freezing (t/f): Previous reported 
values, theoretical prediction from Mean Field Cage Theo- 
ries (MFCT), recent estimations by Van Meel and computed 
values from Polynomial Fitting (PF) and derived Universal 
Relationship (UR). 



D 


Prev. Sim. 


MFCT— 


Van Mee& 


PF 


UR 


3 


0.494±i 


0.494 


0.494 


0.488(5) 


0.492 


4 


0.308 3 


0.308 


0.288 


0.304(1) 


0.308 


5 


0.194^ 


0.169 


0.174 


0.190(1) 


0.189 


6 




0.084 


0.105 


0.114(2) 


0.114 


7 




0.039 




0.0702(2) 


0.0696 


8 




0.017 






0.0427 



expressions to the data in Table IVl are: 

b(D) = 8.942 + 0.00461e 1420 
c(D) = -2.059e- a485D 
d(D) = -0.000945 - 7.082e" 



It is remarkable to note that as D increases, the parame- 
ter b increases fast, indicating that the curve correspond- 
ing to the metastable crystalline states decreases faster 
for higher dimensions and at the limit D — > oo the curve 
changes in a discontinuous way indicating the colapse of 
the freezing and melting point to the same value. 

Equations ((9]) and (fTQ| can also be used to es- 
timate the freezing point, by solving the equation 
gMinhij]-,D) = gMinciv, D). Clearly the solution can 
not be algebraic but anyway provides a semiempirical 
method to determine the freezing point for arbitrary 
dimensions. A list of estimated values for the freezing 
point computed whit the simple the polynomial fit used 
before (PF) and the universal relations ( UR) given by 
Eqs. © and dTUJ) are presented in Table ED 

The results summarized in Table IVII show that the 
use of semiempirical equations improve the estimation of 
the freezing point, at least for D = 4 and D — 5 where 
previous simulations and theoretical approaches are 
available. For D > 6 our estimation always lies above 
the theoretical predictions but no other simulations 
results are available. 

In addition, with the estimated freezing points we can 
extrapolate the melting point and transition pressure us- 
ing the equations of state given in Eqs. ^ and ([7]) for 
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0.0 0.1 0.2 0.3 0.4 0.5 0.6 



n 

FIG. 7. Comparison of numerical simulation values of 
g(Xmin) (dots) with proposed expressions © and (|10[) for 
D=3,...,7. 



TABLE VII. Coexistence reduced pressure (pp) computed 
from the CS equation of state for D = 3 and Pades form 
Bishop et al. for D = 4, 5, 6, 7 



D 


Prev. Sim. 


P*F 

PF 


UR 


3 


11.779 


11.202 


11.668 


4 


11.418 


11.008 


11.469 


5 


14.315 


13.433 


13.184 


6 




16.668 


17.0318 


7 




22.597 


22.1569 



fluid and crystal phases, and using the condition of equal 
pressure in the coexistence region. The resulting values 
are in Tables I VIII and IVIII1 where, although being ex- 
trapolations, may give very good results compared with 
previous simulations up to D = 6. Again an improvement 
in the approach may be remarked from the semiempirical 



TABLE VIII. Melting (t/m) packing fraction (By extrapola- 
tion of p* to crystalline branch) 



D 


Sim. 


t]m 
Van Mee&- 


PF 


UR 


3 


0.545 


0.545 


0.537 


0.542 


4 


0.337 


0.337 


0.368 


0.374 


5 


0.206 


0.206 


0.242 


0.240 


6 


0.138 


0.138 


0.146 


0.147 


7 






0.086 


0.085 



solution. This means that the scheme of chosen equations 
of state and the implicit equation for the freezing point 
is consistent with present available data. 

It is important to remark that the freezing pack- 
ing fractions, estimated using any of the strategies de- 
scribed, change as a function of D following an exponen- 
tial decay, and fit very precisely the function rjp{D) = 
2.103e~°- 482 - D . It is not clear whether it is possible to en- 
sure that this relation will be the same for D out of the 
examined range. In particular for D = 1 and D = 2 
it extrapolates to r] F (l) = 1.298 and rj F (2) = 0.802, 
which do not correspond to the known values (rjp = 1 
for D — 1 and r)p = 0.69 for D = 2). On the other hand, 
the intersection of relations (j9)) and (fT0|) gives the values 
r)F = 1.215 and T]f = 0.756 respectively. Also for D > 7 
the tendency can not be settled as the definitive one. 

At this point, we have presented evidence that the use 
of the properties of the radial distribution function at the 
first minimum can be regarded as an empirical method 
to estimate the freezing point in HHS fluids. However, 
it is not clear why the position and the height of the 
first minimum of the RDF in the crystal and the fluid 
branches could be connected with the freezing point for 
all dimensionalities. The answer is not at all simple, but 
some properties of the three dimensional fluid may give 
some insights. 

In a hard sphere system the fluid-solid transition hap- 
pens, either when the rise or decrease of the free volume 
allows the formation or destruction of cages. The prop- 
erties of such cage formation determine whether the fluid 
form crystals or glasses. In fact, the concept of dynamical 
formation of cages has been already used to approach the 
phase diagrams of hard disks and spheres^. In our case 
our MD simulations began from spheres placed in an FCC 
lattice (for the crystal branch) and high packed random 
states (metastable branch), and we control the packing 
fraction by changing the diameter of the spheres. As is 
usual, this kind of simulations do not allow to achieve 
the coexistence branch because of the finite size of the 
samples^, therefore the crystal branch obtained includes 
a metastable region between the melting and the freez- 
ing points where cages must be present preserving some 
of the crystalline order observed in the RDF. When the 
density decreases the state where those crystalline struc- 
tures melt could correspond to a state where particles 
leave cages, and this happens very close to the freezing 
density. 

The value of X m i n has been usually related with the 
mean number of first neighbors. The integral of the RDF 
in a spherical volume of radius X m i n times the density 
should, in principle, approach this number. In fact, the 
more packaged the system is, the better approach we get 
(at least for the crystal branch). In this context X min 
may be interpreted as the mean radius of a spherical 
volume that contains the first neighbors around a given 
sphere. 

What we observe in Fig. [3] is that the states in the 
crystalline metastable branch approaching the freezing 
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point are still FCC structures where first neighbors are 
confined in average in a sphere of radius lower than y2 
times the diameter, a volume in which an ideal close- 
packed crystal may contain up to the second neighbors. 
It seems that once this volume surpasses this value the 
metastable states quickly "melt". Interestingly enough 
the packing fraction at which this "melting" occurs coin- 
cides with that of the freezing point. 

Recently Kumar and Kumaran^ found, using a 
Voronoi neighbor statistics, the change on the first neigh- 
bor number C\ as a function of the packing fraction for 
the hard-sphere system. They generated hard sphere 
structures using the NVE Monte Carlo (MC) and an- 
nealed Monte Carlo (AMC) methods, to obtain configu- 
rations in two branches, one where crystallization is al- 
lowed and another where not. A plot of C\ as a function 
of r], given in the last reference, shows a behavior similar 
to our Fig. [31 where for the fluid phase C\ has a smooth 
decay when the packing fraction increases, starting from 
values above 15.5 at low densities. Once the number of 
neighbors is around 14.5 the structures obtained by MC 
begin the crystallization and show a sudden decrease in- 
dicating the freezing transition. They obtain a value of 
G\ for crystalline configurations close to 14 instead of the 
expected 12 because a slight perturbation of the position 
in the centers of a pair of molecules may transform a 
vertex in a surface introducing second neighbors to the 
Voronoi count. Annealed structures do not crystallize 
and lead to random dense structures that may be com- 
pared with our metastable fluid branch. 

To inspect the relation among C\ in reference^ and 
X m in we produced Fig. [8] and Fig. [9l relating both vari- 
ables. In Fig. [5] we put together both values assuming 
that the AMC simulations generate configurations in the 
same metastable fluid branch our work does. From the 
figure a proportional relation between both quantities in 
the region of packing fractions close to the freezing point 
is clear. For configurations obtained from NVE (MC) 
one may obtain Fig. [^comparing with the crystal branch 
of our X m i n results, where we note that for comparison 
we moved each function to center both in their respec- 
tive freezing packing fraction, which are close, but are 
different estimations. Despite the differences on simula- 
tion methods it is possible to presume that in both cases 
(more clearly in the metastable branch) there is a linear 
proportion between C\ and X m i n in the vicinity of the 
freezing point. The sudden change on X m i n is propor- 
tional to the change on the number of first neighbors, 
indicating that any remaining cages in the metastable 
states are opening. 

If we put together the last elements we can inter- 
pret that the changes in X m i n and g(X m i n ) close to the 
freezing transition are directly related with the critical 
changes on the cages present in the system. A critical 
number of first neighbors and a critical radius define ei- 
ther the formation or destruction of the cages. In other 
words the dynamical formation of cages, in the context 
cited by Kreamer and Naumis^i, may give the effect ob- 
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FIG. 8. Comparison between values of the first neighbor num- 
ber obtained in^ and the radius of first minimum for the 
liquid-liquid metastable branch. 











1.50 








f 1.45 








* E 1.40 








1.35 


- O 

- o 

t9 







14.1 14.2 14.3 14.4 14.5 14.6 
CArj-rip] 



FIG. 9. Comparison between values of the first neighbor num- 
ber obtained in-^ and the radius of first minimum for the 
liquid-crystal branch. 



served in the first minimum of the RDF. 

For higher dimensions something similar may be hap- 
pening but without a doubt in a more complicated con- 
text. For example, while the D4 and D$ lattices seems to 
be the densest packings for D = 4 and D = 5, in dimen- 
sion D = 6 the Eq lattice is more compact than Dq^L. In 
2006 Skoge et ali21, found that for random jammed states 
the first minimum of the RDF for D = 3, 4, 5 is close 
to where the cumulative coordination number equals the 
kissing number of the densest lattice, giving some clues 
that packing phenomena may present analog properties 
at least up to D = 5. However, in a more recent work 
van Meel et al^ studied the effect of frustration on ran- 
dom packed structures in D = 4 finding that in this sys- 
tem crystallization is slower than in D = 3 because the 
geometrical frustration is surprisingly stronger than in 
D = 3 and that the similarity between the kissing num- 
ber and the number of first neighbors can be explained 
by a wide first peak of the RDF, that accommodates 
nonkissing neighbors in polytetrahedral clusters. Beyond 
these achievements the problem of dynamical formation 
of cages in the freezing and melting transition is not yet 
examined deeply and neither the relation with the first 
minimum of the RDF. In our point of view it is an still 
open problem that deserves some future work. 
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FIG. 10. Comparison of hard-spheres excess entropy as com- 
puted from Eq. (|11[) using values of X m i n obtained from sim- 
ulations (labeled Sim), results from Voronoi cell distribution 
from Ref. |5fJ (labeled Kumar). 



as a consequence, the changes at freezing density are not 
continuous because an entropy gap appears between the 
last metastable crystalline state and the first fluid state. 
Kumar and Kumaran^ report for the entropy gap an 
estimation of As/ Kb ~ 0.92 which corresponds, using 
equation (fTTj) . in the X m i n space to a gap AX TO ,* n ~ 0.15. 
The gap we observe in our simulations is lower than 
Al m j n ~ 0.05. Both of them are small in the scale 
of X m i n and of the order of the error bars, allowing a 
good accuracy for the estimation of the freezing density 
if a continuous change is assumed at the freezing point. 
Since Kumar and Kumaran use the same form of Eq. [TJJ 
for D = 2 and D = 3 with A = D, we believe the same 
situation will apply for D > 3 and Eq. (fTTj) could be 
generalized for arbitrary D (with Vq a function of D) as 

s E (D) 

— M =D\og(V mm (l - V)/V (D)). (12) 
Kb 

This could be the reason why the estimation method for 
the freezing point through X m i n and g(X m i n ) seems to 
work also for the examined dimensions. We plan to per- 
form a further and deeper analysis of the entropy prop- 
erties of hypersphere systems in future work. 



V. SUMMARY AND CLOSING REMARKS. 



In addition it is also possible to trace a path to a ther- 
modynamic interpretation of the observed effects. This 
may be achieved through the examination of the entropy. 
In another recent work of Kumar and Kumaran^, the 
configurational entropy of hard spheres was examined us- 
ing the volume distribution of Voronoi cells. There they 
found a relation between the information entropy of the 
distribution and the excess entropy of the hard-sphere 
fluid, which we have approximately fitted in Fig. 1101 One 
may note that such results can be very well reproduced 
from our results if we define an excess free- volume-like 
entropy per particle as: 

^=AlogflWl-»7)/Vo), ( u ) 

where V rnin = (4Tr/3)X^ lin is the spherical volume corre- 
sponding to X m i n and therefore V m i n (l — if) is the mean 
available volume of a sphere within the first neighbor cell. 
The parameter A is chosen to be A = 3 in accordance 
with the proportionality constant between the informa- 
tion an excess entropies found by Kumar and Kumaran 
and Vq left as free parameter with a value Vq = 40.45. 
From the comparison, shown also in Fig. 1 1 01 it is clear 
that both results must be related and therefore this def- 
inition is fully compatible with the information entropy 
derived from Voronoi statistics. From here we can de- 
rive two conclusions for the system D = 3: first, Fig. 
[3] and Fig. 0] may be taken as two expressions of the 
freezing transition in the entropy-density plane; second, 



We have studied the problem of freezing in HHS flu- 
ids through molecular dynamics simulation from D = 3 
to D = 7, obtaining data for the compressibility factor 
and the RDF in dense liquid states close to freezing and 
crystal metastable states, below the melting density. We 
found that the changes with density of the first mini- 
mum of the RDF, in all examined dimensions, follow a 
tendency where the three branches seems to converge to 
a single point. The packing fraction where this transition 
takes place is a good approximation to the previously es- 
timated freezing points at D = 3, 4 and 5, giving more 
elements to affirm such values are close to the correct 
ones. Assuming that the same structural changes occur 
in higher dimensions we applied the method in our re- 
sults for D — 6 and D = 7, obtaining new estimations 
of the freezing packing fractions that could be in the fu- 
ture confirmed with different theoretical and numerical 
techniques. The proposed method was also formulated 
in a semiempirical scheme, using the RFA method of 
Rohrman and Santos^ to obtain the first minimum of 
the RDF numerically in the fluid phase and finding the 
intersection with empirical relations in the crystal phase. 
General relations for the value of the RDF in the first 
minimum, as a function of the dimension and the pack- 
ing fraction were proposed, giving also very good approx- 
imations to known values for the freezing points. 

Although the sets of particles used in this work are 
small, we are confident that the results describe the over- 
all system behavior. Previously it was noted that for di- 
mension D = 7, the equation of state is well described 
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even for a small set of particles (N = 64) . A simple com- 
parison of the g(r) properties shows that the influence 
of the system size is reflected as a slower convergence 
and a more noisy numerical result, but the mean value of 
g{X m in) seems to be consistent for different sizes of the 
system. However, a more detailed analysis may clarify 
this conjecture. 

An entropy analysis reveled that the properties of 
X m i n and g(X m i n ) may be related with the changes in 
the excess entropy of the system and that a small gap be- 
tween the fluid and the crystalline metastable branches 
exists, but the size of it is small enough to allow the use of 
the proposed numerical method to obtain the estimations 
of the freezing density with good accuracy. 

Finally, it is worth to mention that the method of 
tracking the minimum of the RDF is reliable for the ex- 



amined cases, but as we increase the dimensionality there 
is an effect of flatness in the fluid RDF and sharpness in 
the the crystal RDF, collapsing all minima close to 1 in 
the fluid phase and close to in the crystal phase, and 
therefore, the accurate determination of the position and 
height of either the maximum or minimum becomes more 
difficult. 



ACKNOWLEDGMENTS 

We acknowledge the financial support of DGAPA- 
UNAM through the project IN-109408-2. One of us, C. 
D. E. acknowledge to CONACyT for financial support 
through the grant 171940. We also want to thank Prof. 
Mariano Lopez de Haro for helpful comments. 



* [cdeaQcie.unamjmx] 
t mrp@cie.unam.mx 

1 F. H. Ree and W. G. Hoover, J. Chem. Phys. 40, 2048 
(1964). 

2 E. Leutheusser, Physica A 127, 667 (1984). 

3 J. P. J. Micfiels and N. J. Trappeniers, Physics Letters 104, 
425 (1984). 

4 B. C. Freasier and D. J. Isbister, Mol. Phys. 42, 927 (1981). 

5 M. Luban and A Baram, J. Chem. Phys. 76, 3233 (1982). 

6 C. G. Joslin, J. Chem. Phys. 77, 2701 (1982). 

7 Theory and simulation of Hard Sphere fluids and related 
systems, edited by A. Mulero, Lecture Notes in Physics, 
Vol. 753 (Springer, 2008). 

8 R. Finken, M. Schmidt, and H. Lowen, Phys. Rev. E 65, 
016108 (2001). 

9 B. J. Alder and T. E. Wainwright, J. Chem. Phys. 27, 1208 
(1957). 

10 B. J. Alder and T. E. Wainwright, J. Chem. Phys. 33, 1439 
(1960). 

11 B. J. Alder and T. E. Wainwright, J. Chem. Phys. 31, 459 
(1959). 

12 W. W. Wood and J. D. Jacobson, J. Chem. Phys. 27, 1207 
(1957). 

13 W. G. Hoover and F. H. Ree, J. Chem. Phys. 49, 3609 
(1968). 

14 E. G. Noya, C. Vega, and E. de Miguel, J. Chem. Phys. 
128, 154507 (2008). 

15 Xian-Zhi Wang, J. Chem. Phys. 122, 044515 (2005). 

16 P. Tarazona, Phys. Rev. A 31, 2672 (1985). 

17 W. A. Curtin and N. W. Ashcroft, Phys. Rev. A 32, 2909 
(1985). 

18 A. R. Denton and N. W. Ashcroft, Phys. Rev. A 39, 4701 
(1989). 

19 J. F. Lutsko and M. Baus, Phys. Rev. A 41, 6647 (1990). 

20 Yaakov Rosenfeld, Phys. Rev. Lett. 63, 980 (1989). 

21 Z. Cheng, P. M. Chaikin, W. B. Russel, W. V. Meyer, J. 
Zhu, R. B. Rogers, and R. H. Ottewil, Materials & Design 
22, 529 (2001). 

22 G. Bryant, S. R. Williams, and L. Qian, Phys. Rev. E 66, 
060501 (2002). 

23 M. Robles, M. Lopez de Haro, and A. Santos, J. Chem. 
Phys 120, 9113 (2004). 



24 D. J. Gonzalez, L. E. Gonzalez, and M. Silbert, Mol. Phys. 
74, 613 (1991). 

25 M. Bishop and P. A. Whitlock, J. Stat. Phys. 126, 299 
(2007). 

26 M. Bishop, N. Clisby, and P. A. Withlock, J. Chem. Phys. 
128, 034506 (2008). 

27 M. Skoge, A. Donev, F. H. Stillinger, and S. Torquato, 
Phys. Rev. E 74, 041127 (2006). 

28 J. A. van Meel, D. Frenkel, and P. Charbonneau, Phys. 
Rev. E 79, 030201 (2009). 

29 J. A. van Meel, B. Charbonneau, A. Fortini, and P. Char- 
bonneau, Phys. Rev. E 80, 061110 (2009). 

30 T. M. Truskett, S. Torquato, S. Sastry, P. G. Debenedetti, 
and F. H. Stillinger, Phys. Rev. E 58, 3083 (1998). 

31 J. P. Hansen and L. Verlet, Phys. Rev. 184, 151 (1969). 

32 H. R. Wendt and Farid F. Abraham, Phys. Rev. Lett. 41, 
1244 (1978). 

33 M. Bishop, P. A. Whitlock, and D. Klein, J. Chem. Phys. 
122, 074508 (2005). 

34 C. D. Estrada, Study of the system of Brownian hard 
spheres,, M. sc. thesis, Posgrado en Ciencias Ffsicas, Uni- 
versidad Nacional Autonoma de Mexico (2005). 

35 N. F. Carnahan and K. E. Starling, J. Chem. Phys. 51, 
635 (1969). 

36 K. R. Hall, J. Chem. Phys. 57, 2252 (1972). 

37 R. J. Speedy, J. Chem. Phys. 100, 6684 (1994). 

38 M. N. Bannerman, L. Lue, and L. V. Woodcock, J. chem. 
Phys. 132, 084507 (2010). 

39 S. B. Yuste and A. Santos, Phys. Rev. A 43, 5418 (1991). 

40 A. Trokhymchuk, I. Nezbeda, J. Jirsak, and D. Henderson, 
J. Chem. Phys. 123, 024501 (2005). 

41 L. Lue, M. Bishop, and P. A. Whitlock, J. Chem. Phys. 
132, 104509 (2010). 

42 M. Bishop and P. A. Whitlock, J. Chem. Phys. 123, 014507 
(2005). 

43 G. Nebe and N. J. A. Sloane, "A catalogue of lattices," 
http://www2.research.att.com/~njas/lattices/. 

44 R. J. Speedy, J. Phys. C 10, 43874 (1998). 

45 M. Luban and J. P. Michels, Phys. Rev. A 41, 6796 (1990). 

46 R. D. Rohrmann and A. Santos, Phys. Rev. E 76, 051202 
(2007). 



13 



A. S. Kraemer and G. G. Naumis., J. Chem. Phys. 128, 
134516 (2008). 

V. S. Kumar and V. Kumaran, J. Chem. Phys. 123, 074502 
(2005). 



J. H. Conway and N. J. A. Sloane, Sphere packings, lattices 

and groups (Springier- Verlag, New York, 1998). 

V. S. Kumar and V. Kumaran, J. chem. Phys. 123, 114501 

(2005). 



